home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QLOG.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  1KB  |  87 lines
  1. /*                        qlog.c    */
  2. /* natural logarithm */
  3. #include <stdio.h>
  4. #include "qhead.h"
  5.  
  6. extern QELT qone[], qtwo[], qlog2[], qsqrt2[];
  7.  
  8. int qlog( x, y )
  9. QELT *x, *y;
  10. {
  11. QELT xx[NQ], z[NQ], a[NQ], b[NQ], t[NQ], qj[NQ];
  12. long ex;
  13.  
  14. if( x[0] != 0 )
  15.     {
  16.     qclear(y);
  17.     printf( "qlog domain error\n" );
  18.     return 0;
  19.     }
  20. if( x[1] == 0 )
  21.     {
  22.     qinfin( y );
  23.     y[0] = -1;
  24.     printf( "qlog singularity\n" );
  25.     return 0;
  26.     }
  27. /* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
  28. qmov(x, xx );
  29. ex = *(xx+1);
  30. if( ex == EXPONE )
  31.     { /* log 1 = 0 */
  32.     if( qcmp(x, qone) == 0 )
  33.         {
  34.         qclear(y);
  35.         return 0;
  36.         }
  37.     }
  38. ex -= (EXPONE-1);
  39. xx[1] = (EXPONE-1);
  40. /* Adjust range to 1/sqrt(2), sqrt(2) */
  41. qsqrt2[1] -= 1;
  42. if( qcmp( xx, qsqrt2 ) < 0 )
  43.     {
  44.     ex -= 1;
  45.     xx[1] += 1;
  46.     }
  47. qsqrt2[1] += 1;
  48.  
  49. qadd( qone, xx, b );
  50. qsub( qone, xx, a );
  51. if( a[1] == 0 )
  52.     {
  53.     qclear(y);
  54.     goto bdone;
  55.     }
  56. qdiv( b, a, y );    /* store (x-1)/(x+1) in y */
  57.  
  58. qmul( y, y, z );
  59.  
  60. qmov( qone, a );
  61. qmov( qone, b );
  62. qmov( qone, qj );
  63. do
  64.     {
  65.     qadd( qtwo, qj, qj );    /* 2 * i + 1        */
  66.     qmul( z, a, a );
  67.     qdiv( qj, a, t );
  68.     qadd( t, b, b );
  69.     }
  70. while( ((int) b[1] - (int) t[1]) < NBITS );
  71.  
  72.  
  73. qmul( b, y, y );
  74. y[1] += 1;
  75.  
  76. bdone:
  77.  
  78. /* now add log of 2**ex */
  79. if( ex != 0 )
  80.     {
  81.     ltoq( &ex, b );
  82.     qmul( qlog2, b, b );
  83.     qadd( b, y, y );
  84.     }
  85. return 0;
  86. }
  87.